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ABSTRACT 

We present spline-based techniques for estimating spatially varying 
parameters that appear in parabolic distributed systems (typical of those 
found in reservoir simulation problems). In particular, we discuss the 
problem of determining discontinuous coefficients, estimating both the 
functional shape and points of discontinuity for such parameters. In 
addition, our ideas may also be applied to problems with unknown initial 
conditions and unknown parameters appearing in terms representing external 
forces. Convergence results and a summary of numerical performance of the 
resulting algorithms are given. 
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1. Introduction 


We present here our efforts related to the estimation of discontinuous 
spatially varying coefficients in parabolic distributed systems. Although our 
ideas are applicable to a wide class of problems in which the determination of 
discontinuous coefficients is of importance (e.g., the propagation of waves 
through layered media; the dynamics of beams with "discontinuous" elastic 
properties), our work here is motivated by an inverse problem in reservoir 
simulation commonly referred to as "history matching". The problem in this 
case is to determine unknown parameters (such as permeability, porosity) that 
appear as coefficients in model reservoir equations. "Optimal" choices of these 
parameters should provide the best match between the observed and simulated 
production history at one or more wells. Information about these coefficients 
(functional shape and location of discontinuities) provides insight into physical 
properties of the reservoir and can indicate the location of abrupt structural 
changes; in addition, precise determination of these parameters is essential to 
the process of accurately simulating and predicting reservoir behavior. 

The governing reservoir equations describe mathematically the physical 
and chemical processes occurring during primary hydrocarbon recovery or during 
enhanced recovery efforts (secondary or tertiary forms of recovery). Mathematical 
models vary widely depending on the physical process being described (miscible 
or immiscible fluid flow, thermal or fluid injection, etc.) and the types 
of observations available. Common to each model however is a system of rate 
equations (derived from Darcy’s law, which relates flow rate to fluid pressure 
gradients) as well as appropriate conservation laws and equations of state. 

The resulting dynamical system is typically distributed in nature and of para- 
bolic type [17], [18]; unknown parameters quite often include the porosity of 
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surrounding rock, or the ratio of pore volume to total volume, and (relative) 
permeability, which is the ability of the rock to transmit fluid [18]. Due to 
spatial changes in underground structure, it is highly likely that these 
parameters will vary spatially and contain numerous discontinuities. 

In order to solve the inverse problem, data in the form of fluid pressure 
(or flow rate) is collected at the wells and used in a numerical parameter 
estimation process. There have been a large number of substantial contributors 
to the development of theoretical concepts and numerical algorithms for the 
history matching problem. An exhaustive list of related references would be 
too lengthy to include here; instead we refer the reader to [18] for an excellent 
survey of the outstanding efforts in this area. One numerical approach commonly 
taken involves subdividing the reservoir into a grid of smaller blocks; constant- 
valued parameters (which are allowed to vary independently from block to block) 
are then estimated. Unfortunately, if accurate solutions are desired, the 
grid size often must be quite small and thus the number of unknown parameters, 
as well as the dimension of the state space, can be very large--as many as 
50,000 parameters or more [17]. (This is an unfortunate consequence of the 
fact that the parameters of interest — as well as state variables — are infinite- 
dimensional yet computations must be performed in a finite-dimensional setting.) 
Our goal here is to avoid some of the difficulties associated with the approach 
described above. Specifically, our ideas involve separating the order of 
state approximation from that of parameter estimation, so that the need for 
an approximate state space of high dimension does not impose the same require- 
ments on the dimension of an approximate parameter space; this is accomplished 
by searching for parameters in classes of functions with quite general spatially 
varying representations. In order to focus attention on the problems associated 
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with estimating spatially varying discontinuous coefficients in this context, 
we consider an archetypical model of (parabolic) distributed type that 
admittedly is a simplified version of the fluid pressure equations associated 
with reservoir simulation (see [17], [18], and the references therein); never- 
theless the model selected here is a prototype that contains the essential 
parameter-dependent terms for which we may begin our investigations. In the 
sections that follow we define the model equations of interest and construct 
an approximation framework in which we wish to consider the parameter estimation 
problem. Convergence results are presented for problems associated with either 
spatially distributed or "discrete" sample data. Finally, we discuss numerical 
implementation in general, and in the context of particular examples. 


It is our intent in this report to examine convergence properties and imple- 
mentation problems associated with these methods; we do not address such important 
questions as identifiabili ty, observability, or general underlying properties 
of the governing partial differential equation system. 


The notation used throughout is standard: For I C R (the real line), we 

shall denote by C ( I ; X ) the space of continuous functions f : I -*■ X with uniform 
norm | - 1 ^ ; by L 2 (I;X) we mean the usual space of square-integrable "functions" 


and inner product <*,*> L (i-x)‘ The Sobolev 


f : I -> X wi th L 2 norm | • | L ( j . x ^ 

spaces H^(I;X) and Hq(I;X) are defined as usual (see, for example, [i ]). 

Whenever X = R, we shall simplify notation by writing C(I) and l_ 2 (I), respectively, 
and, where no confusion results, by writing |*| (and <*,•>) for the norm (and 
inner product) on L 2 (0,1). In addition, no notational distinction will be 
made between a function f : I -> R and its restriction to ^ C I. 
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2. The Parameter Estimation Problem 

As our fundamental state system we consider the scalar parabolic distributed 
system 

^ (t - x) « (q(x) If (t,x)) + f < t - x ; r < x >>- (°> T > * (0 ' 1) - 

5 ] ^ u(t,0) = u(t,l ) = 0 , 

u(0,x) = Uq(x) . 

Here q and p are discontinuous (positive) functions representing the permeability 
and porosity properties, respectively, of the fluid and surrounding rock; the 
points of discontinuity in these functions correspond to abrupt spatial changes 
in the physical flow region (such as might be associated with layered media). 

Both q and p are typically unknown so we shall consider the problem of estimating 
these parameters, as well as the function r, r(x) e R and the initial condition 
Uq, from observations of the state variable u. 

For ease of presentation in the arguments that follow, it is assumed that 
p = 1, although it is not difficult to extend our ideas to the case of non- 
constant (and unknown) p. We detail in Remark 3.1 the minor modifications one 
must make in the calculations found below in order to treat p as a functional 
parameter throughout. 

To simplify notation, we assume that q is discontinuous at one point only, 
x = ?, and that q is represented by 


q = 
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where <j>-j and cf> 2 are continuous on [0,1]; here is the usual Heaviside function 
on [0,1] given by H^ = 1 on [5,1], H = 0 otherwise. There is a straightforward 
extension of our ideas to the case where 
y 

<1 = 'f’l + I H r <f>- 

1 i=2 5 i-l 1 

0 = 5q < 5-j < 5 2 < • • • < - 1 , except that notational difficulties become 

excessive. (We later demonstrate our approximation and estimation techniques 
for multiple discontinuity problems in the section on numerical findings.) Given 
the parameterization chosen for q we define the parameter vector y = (5, <t>-| > <f>2» 
r, u q ) = (s,Uq) as an element of the parameter set r c r = 3 x Lg (0 , 1 ) , where, 
for m, m fixed, 

S(m,m) e {s = (Z,<p r <p 2 ,r)£ R xC[0,l] x C[Q,1 ]. * L 2 ( (0,1) ;R P ) | 5 .-1(0,1), 
<j>. £ ^[0,1], and 0 < m £ <f>.j(x) £ m 
for i = 1,2, and x 6 [0,1 ]}. 

Concerning r and the applied force f, we make the following (standing) hypotheses: 
(HI) The parameter set r is compact; 

(H2 ) For every r € L 2 ((0,1 );R P ), the map t -* f(t,-;r(-)) : [0,T] + L 2 (0, 1 ) 
is Holder continuous with exponent a, 0 < ct < 1. 

(H3) The map r f(* ,* ; r ( - )) is continuous from L ((0,1); R P ) to L 2 ((0,T) 

(0,1)). 


X 
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The parameter estimation problem associated with (2.1) consists of finding a 
★ 

parameter y e r that is "optimal" in the sense of providing the best match 

between observed data and model solutions to (2.1). Although a number of 

criteria may be used to measure "fit to data," we consider first a least squares 

criterion J that is defined in conjunction with distributed data: That is, 

given distributed observations u^ G l^Ojl) at discrete times t^ € (0,T), 

★ 

i = l,...,n, we seek y e r that minimizes 

n 1 o 

(2.2) J(y) = l / |C(t.,x;y)u(t.,x;y) - u.(x)| dx 

i=l 0 1 1 1 

over all y G r. For each (t. ,x), the output map C(t. ,x;y) : R R is assumed 

to be continuous in y and such that the mapping x -* C(t^ ,x;y)ij;(x) is in l^O.l) 

whenever \p ^ Lg (0,1 ) . We note that data generally is not available in the 
distributed form given here; often this difficulty can be handled by fitting 
a curve (using linear interpolation, for example) to discrete data. 

We also treat the problem of truly discrete data, i.e. , u..G R is observed 

* J 

sample data at (t.,x.), j = l,...,n. In this case the parameter estimation 

J 

~ ★ 

problem consists of determining y € r that minimizes a "pointwise" fit- to- data 
criterion, 

(2.3) J(y) = l l |C(t.,x.;y)u(t.,x.;y) - u..| 2 

i=l j=l 1 J 1 J 
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over y 6 r. The use of discrete sample data leads to increased technical detail 
and additional smoothness hypotheses on Uq and f. We defer consideration of 
this particular estimation problem until we have fully developed an approximation 
theory for distributed estimation problems (i.e., identification problems with 
cost functional J defined in (2.2)); our findings for the "pointwise" estimation 
problem (involving J) are then summarized in section 3.1 below. 

Before we consider either of these estimation problems (where y£ r is 
unknown and to be determined), we shall first consider the existence of solutions 
u of (2.1) for a given parameter y = ( 5 , 4 >^ , 4 > 2 »r,UQ) ^ r. Defining 
u(t) = u(t,*) € LgtO.l), we may rewrite (2.1) as an initial value problem in u, 


(2.4) 



= A(q)u(t) + F(t;r) , 




U(0) = Uq . 


t € (0,T) , 


Here q = + H F(t;r) = f(t,*;r(*)) and the operator A(q) is defined by 

A(q )ip = D(qDi|») for ip € domA(q) = V , where 

v q = 4 e Hq(0, 1 ) I qDip e H 1 (0,1 )} 

(throughout we shall use D to denote the spatial differentiation operator — ). 
We note that it would be natural, given the discontinuities in coefficients p 
and q, to consider a weak form of (2.1) in order to relax restrictions on both 
solutions and parameters. For this particular problem, however, we shall 
insist that solutions u satisfy a continuity equation 


(2.5) (qDu)U') = (qDu)U + ) , 


which represents continuity of stress across a transition point, £, between 
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distinct spatial regions (layers of porous media, for example). Given this 
condition on u (which implies that qDu must be sufficiently regular to ensure 
that point evaluations are meaningful), there is reason for seeking at least 
strong solutions to (2.1). We have one additional comment along these lines 
from a numerical point of view: It is of interest to note that we experienced 
no difficulties in applying our state variable approximation ideas to the 
forward problem (i.e., the problem of integrating (2.1), or a weak form of 
(2.1), for a known value of y) for particular problems where the true state u 
did not satisfy (2.5). It was not until we turned to the numerical solution 
of the i nverse problem (and, in particular, the problem of estimating £ itself) 
that it became evident that one could not expect to estimate £ unless the 
(physically meaningful) continuity equation (2.5) was satisfied by solutions 
of (2.1). Therefore, there seems to be little justification in considering weak 
solutions of (2.1) in the context of the parameter estimation problem. 

Our first result is a statement of existence and uniqueness of solutions 
of (2.1); in addition, we indicate certain regularity properties of solutions 
that will be useful in later calculations. 


Theorem 2.1. Let y = (g,<f>.| ,cj>2,r,Ug) be given in r and let q = ^ + H $ 
There exists a unique (classical) solution u to (2.1) with the property 
that u(t)€ 1/^ for any t > 0. In addition, if Uq € V , then the map 
t -* A(q)u(t) is in C([0,Tj; L 2 (0,1 )). 


2 ’ 


Proof : It suffices to show that A(q) generates an analytic semigroup on 

L 2 ( 0 , 1 ) . From this we may guarantee existence of a unique solution u to (2.1) 

(see Corollary 3.3, p. 113 of [27]); we may then apply the well-known smoothing 
properties of analytic semigroups (see, for example, Theorem 3.5, p. 114, of [27]), 
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along with hypothesis (H2), to obtain the statement of the remainder of the 
theorem. 

We first show that A(q) is densely defined and self-adjoint. To this end 

we note that = {> € L 2 (0,1) | ip G Hq(0,?), Hq(£, 1)} satisfies 

V q C V q = domA(q) (since qD^(^") = q!tyU + ) = 0); it is easy to argue that V q 

2 2 

is dense in l_ 2 (0,l) if one uses the fact that HQ(0,?)and Hq(£, 1) are dense in 
L 2 (0, c ) and L 2 U,1), respectively. Using an integration by parts, it is not 
difficult to show that A(q) is symmetric. To demonstrate that A(q) is self- 
adjoint, it suffices to show [28; Theorem 13.11] that Range A(q) = L 2 (0,1); that 
is, for g € L 2 (0,1), it is sufficient to- verify the existence of a solution to 

(2.6) A(q)^ = g 

that satisfies ^ £. V . Because one may readily see this is true (using standard 
theory for two-point ordinary differential equation boundary value problems — 
see, for example. Theorem 8.3 of [19]), it therefore follows that A(q) is self- 
adjoint. In addition, A(q) is dissipative (since <A(q)<ji,i|>> £ -m | Dip | £ 0 for 

all tJi € V ) so that a corollary of the Lumer Phillips theorem [27; p. 15] may 
be invoked to argue that A(q) generates a Cq semigroup of contractions on 
L 2 (0,1). Finally, since a(A(q)) C (-»,0), we may apply standard semigroup 
theory [20; Theorem 7.12, p. 82], [27; p. 61] to conclude that A(q) generates 
an analytic semigroup (analytic on the sector {\G C | A f 0, |argA | < j }). 

The proof of the theorem is complete. 

It is useful to note that if u is a solution of (2.1) then u also satisfies 
(2.1) in a weak sense; i.e., u satisfies 
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(2.7) 


<u t (t),v> 


u(0) = u 


0 


<qDu(t) ,Dv> + <F(t;r),v>, 


t € (0,T) 


for any v£ Hg(0,l). Using this formulation we may argue the continuous 
dependence of solutions on (possibly unknown) initial data. In fact, we may 
actually show that the map y u(t;y) : r -*■ L2 (0,1) is continuous, uniform in t, 
so that we are guaranteed the existence of a minimizer for J over the (compact) 
set r. As we also establish existence of an "optimal" parameter in Theorem 3.3 
(and because we do not need the continuity of y -»■ u(t;y) to make any of the 
arguments given below), we state and prove only the result that follows. 


Corollary 2.1 . The mapping Ug -> u ( t ; 5 ,<j> 2 *r,Ug) : L 2 (0,1) ■> L 2 (0,1) is 
continuous, uniform in ( c , 4>-j ,<f> 2 ,r) € $ and t 6 (0,T). 

Proof : Let Uq, Uq be given in L 2 (0,1) and define u(t) 5 u(t;?,<j>^ , 4 > 2 ,r,UQ) , 

u'(t) = u(t; 5 ,<j)-j ,(|) 2 ,r,Ug) . Using (2.7), we find that, for t>0, 

<u^(t) - uj.(t),v> = - <qD(u(t) - u'(t)),Dv> 

for any choice of test function v in setting v = u(t) - u'(t)e Hg (from 
Theorem 2.1) we obtain 

" u' (t )| 2 + m | D(u(t) - u‘(t))| 2 < 0 . 


In fact, using the Rayleigh-Ritz inequality [29; p. 5], it follows that 
\ |u(t) - u'(t)| 2 + mu 2 | u( t) - u'(t)| 2 < 0 

so that an application of the Gronwall inequality yields 
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I u { t ) - u'(t)| 2 < e- 2 "" * |u 0 - u^l 2 . 


Continuous dependence of u(t;s,u Q ) on u Q , uniform in ses, thus obtains. 


We are ready to consider the problem of unknown parameters, in the context 
of the parameter estimation problems defined in this section. We note that 
the problem of estimating an optimal parameter y € r must be combined with 
schemes for solving (2.1); i.e. , we must consider state variable approximation 
as well as the problem of finding finite-dimensional approximations for unknown 
functional parameters. In the sections that follow, we develop both state and 
parameter approximation techniques with the goal of solving these problems. 
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3. A Spline-Based Approximation Scheme 

Standard numerical optimization schemes applied to the problem of minimizing 
J (or J) over r typically generate a minimizing sequence of parameter iterates, 
starting from an initial guess, y°. However, schemes of this type generally 
require that u(y) (the solution of (2.1)) be evaluated as the parameter y is 
updated; it is therefore desirable to combine estimation of an optimal parameter 
y* with approximation techniques for solving (2.1). With this goal in mind, 
we describe a spline-based state/parameter approximation scheme in the same 
spirit of the ideas found in [5], [9], [10], [12], [14], [23] to name a few 
of the related references in this area for (continuous coefficient) parabolic 
problems. 

The convergence arguments developed below are similar to standard variational- 
type estimates often used in association with finite element approximations (see, 
for example, [29], p. 129) although the estimates given here are complicated 
somewhat by the presence of unknown parameters. This variational approach was 
taken in [10] and [11] for the problem of estimating continuous coefficients in 
parabolic systems; we require a somewhat different treatment here primarily due 
to the fact that we allow discontinuous coefficients, where the points of dis- 
continuity are unknown (necessitating parameter-dependent approximation spaces 
X^*(q)). Thus, an interesting aspect of our approach (and often a source of 
difficulties) involves the fact that our approximation spaces change with every 
choice of parameter iterate. We note that although the theoretical problems 
are quite different, our construction of approximating spaces X^(q) is somewhat 
similar to the ideas found in [2], [13], [16]; there the problem was to esti- 
mate unknown delays appearing in functional differential equations 
(there is a correspondence between our treatment of an unknown point of dis- 
continuity and the approach taken in those references to handle an unknown 
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delay, at least from the standpoint of numerical approximation schemes). 

We turn now to a precise statement of the approximation scheme under 
consideration. 

For any y = (£,<j>-| r » u 0 ) €. r, we construct parameter-dependent spaces 

and operators as follows: For q = ^ + H^ 2 and N = 1,2,..., we define 

X (q) = spantB^q), i = 1,...,2N-1}, where B^(q) denotes the i^* 1 continuous 
piecewise-linear B-spline basis element (satisfying homogeneous boundary con- 
ditions) with knots at (x^(q), k = 0 2N}. Here xjj(q) = k£/N, k = 0,...,N, 

and x^(q) C + ( k- N) (l-c)/N, for k = N+1,...,2N. The piecewise linear 

elements are characterized by B^(q)(x^) = <5^ for i,k = 1 ,2N-1 

(B^(q)(0) = B^(q)(l) = 0); see Figures 1-3. 



Figure 1. bJ, i = 1,...,N - 1. 



Figure 2. b{J . 



14 



Figure 3. B^, i = N + 1 , . . . , 2N - 1 . 


We remark that in general, for y, y£ r and q = ^ q = d>-| + H~<j> 2 , we 

do not have X^(q) C X^(.q) or X^(q) C X^(q), nor do we have X^(q) C (note 

that although an element ip^E X^(q) does have a discontinuity in its first 
N 

derivative at £, \p does not satisfy the continuity equation (2.5) 
associated with q. The approximation spaces X (q) are chosen so 
that the resulting parameter estimation algorithm enjoys a number of computational 
advantages, especially when 5 is unknown. We take a general Galerkin approach 
to define approximating state systems and then obtain convergence findings by 
working directly with the weak form of these equations. As an alternative 
approach to that taken here one could define approximating operators A^(q) for 
A(q) and investigate the sense in which A^(q) "converges" to A(q) (see, for 
example, Example 2.2 of [22], or [23], for approximating operators that might 
be used in this context). 

For ye r fixed and N = 1,2,..., we seek an approximation to u(t;y) of the 

N 2N-1 N N N 

form u N (t;y) = [ w. (t;y)B. (q) , where the "Fourier coefficients" w. are 

i=l 1 1 1 

determined via the system of ordinary differential equations (ODE), 
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(3.1) 


<uj(t;r), B^(q)> = - <qDu N (t;y), DB^(q)> + <F(t;r), B^(q)>, t£ (0,T), 


<u N (0;y), B^(q)> = <u Q , B^(q )> , 


N 

for i = 1,...,2N-1. Alternatively, u satisfies 
(3.2) 


<u!) ( t,y ) , v> = - <qDu N (t;Y),Dv> + <F(t;r),v> , t € (0,T), 


v u N (0;y) = P N (q)u Q 


for all v€ X N (q); here P^(q) : [^((M) -> X^(q) denotes the orthogonal projection 

(with respect to the usual Lg topology) along X N (q)' ! *. Associated with (3.1) 

-N 

is an approximate estimation problem, namely that of finding y G r that minimizes 


n 1 

(3.3) J N (y) =11 |c(t.,x;y)u N (t.;y)(x) - u,(x)| 2 dx 
i=l 0 1 1 1 

over r, where u^(y) is the solution of (3.1) corresponding to y € r. 

X L. 

Our initial findings concerning the N tn approximate problem (3.1), (3.3) 
are immediate consequences of the fact that (3.1) is an ODE on X N (q) and that 
the basis elements B^(q) (and their spatial derivatives) are continuous in 5 
(see (4.1) for a more explicit matrix representation of (3.1)). We shall defer 
a more detailed examination of (3.1) (and the implementation of the estimation 
scheme associated with (3.1), (3.3)), until Section 4 where our numerical 
findings are summarized. 

Theorem 3.1 . For each N and any y 6 r, there exists a unique solution u^(y) 
of (3.1), u N (t;y) e X N (q). In addition, the mapping y -> u N (t;y) : r -> Lg ( 0 , 1 ) 
is continuous for each t 6 (0,T). 
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-N 

Corollary 3.1 . For each N, there exists a solution y £ r for the problem of 
N 

minimizing J over r. 


Finally, a simple modification of the proof of Corollary 2.1 yields a 
similar statement concerning the continuous dependence of u^(t) on Ug: 

Corollary 3.2 . The mapping Ug u^(t;(s,Ug)) : [^(Ojl) -> [^(Ojl) is con- 
tinuous, uniform in N, s=(£, ,j> , , r ) € 5, and te (0,1). 


An essential step in the process of correlating state variable approximation 

* k 

with the problem of estimating an optimal parameter y 6. r (for the original 

parameter identification problem) is the establishment of the convergence of 

u^(t;y^) to u(t;y) for any sequence {y N } in r that converges to y € r. We shall 

clarify the need for arguments of this type in the proof of Theorem 3.3. To 

facilitate steps in this direction, we shall first establish linear spline 

estimates, the proof of which are in the spirit of [29; pp. 16-17, 78], In 

what follows we assume that {y N } is given in r, y^ = j^.r^.Ug) , with 

y^ 1 -* y = (£,<t>i j^jf'jUg) G r (in the usual product topology on r); in addition, 

we assume there exists 6 such that 0<6<_i;<l-5<l (and, in the case of 

multiple discontinuities, - S k _-| | ^.6 > 0, k = l,...,y). Given 
N N N 

q = <)>.| + H we shall henceforth simplify notation and abbreviate 

D N _ D N, N, 5 y N _ Y N/ N, . N _ N/ N, . _ n 
P - P ( q ) 5 X — X (q ) 5 and x. = x.(q ) 5 k 0 9 *«. 5 2N. 


Lemma 3.1 . Let y be given in V-, where q 
c^ and Cg, independent of N, such that 

(3.4) In. - P N ^| < Cl N" 2 |A(q)^| 


4>1 + H-^- There exist constants 
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and, for N sufficiently large, 

(3.5) |D(* - P N *)| < c 2 N _1 |A(q)#|. 

Proof. We shall denote by 1% the linear interpolant of ^ , with knots at 
x k* k ~ 0»-‘*>2N; that is, I ^(xj^) = k = 0,...,2N. We find that 

|D(^-1%)| 2 = <Xty,p(^-i%)> - <D(l%),D(^-i N ^)> 

where N 

/ N n 2N k M .. 

<D(l ip), D(tfj-I \p)> = 1 f D(1 

k=l N 
x k-l 


2N 

I / 

k=l x 


N 

k-1 


D (I N ^) ) 


= 0 . 


Thus, 


|D(^-I N ^) I 2 = <D^,D(*-1%)> 

1 2N X k - |\i 

1 ml / W 0(^-1 

k=1 x N 
x k-l 

1 2N X k - M 

■“^1 / O(qD^) (ip-x ijj) 

k =l vN , 

x k-l 

|p(w)| |lVt| 

i (rmtN )' 1 |A(q)*)| |»(*-I N *)| 
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where we have used (2.16) of [29; p. 17] in the last estimate. It therefore 
follows that 

(3.6) | D(ip-l N ip) | £ (m-irN ) 1 |A(q)<Ji| 


and, again using (2.16) from [29], 

(3.7) |lV*| 1 (m7r 2 N 2 ) -1 |A(q)*| . 

N 

To establish (3.4), we use properties of the projection P to note that 
k- p %l 1 1 c^N 2 |A(q)i|>| • 

Finally, |D(i|i-P%)| £ |d(^-X N iJj)| + |D(I%-P%) | , where an application of the 
Schmidt inequality [29; p. 7] yields (for N sufficiently large) 

HXxVpW- T /*" IiXiVpW 

k=l X N . 


2N 

12 I 

k=l 


(x^-x N 


. -2 


k-1 


) / I 

x N 

X k-1 


x\ 


■p%l 2 


< 12 (N/6) 2 |lVp N t|/| 2 


so that 

j D(ip-P N ^) I £ | D(ip-I N ^) | + 2/3~(N/6) |l%-P%| 

£ | D ( ip— ) | + 2/3" (N/6 ) { | | + |i|j-P%|} 

£ (iattN) ^ | A( q ) ip | + 4/3"(N/5)c-|N 2 [ A(q )^ j j 

£ c 2 N _1 | A(q )ip | 


for an appropriate choice of the constant c 2> 
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We may now use the linear spline estimates derived in Lemma 3.1 to 
establish a preliminary convergence result. 

Lemma 3.2 . Suppose {y N } is given in r such that y N y = (i,^ ,$ 2 ,r,u Q ) £ r. 
Assume, in addition, that Uq £ V- , where q = ^ + H|<j> 2 . Then, for every 

t e (o,t) , 

u N (t;y N ) u(t;y) in L 2 (0, 1 ) 

as N + * (where u N is the solution of (3.1) associated with y N and u is the 
solution of (2.1) associated with y). 

Proof. Let u(t) = u(t;y), u N (t) = u N (t;y N ). Then 

|u N (t) - u ( t ) | < |u N (t) - P N u(t) | + |P N u(t) - u ( t) | 

where the second term is 0(N -2 ) from (3.4) (u(t)€ V- for t£ (0,T); Theorem 2.1) 
To consider the first term, we note that solutions u(t), u N (t) of (2.1), (3.1), 
respectively, satisfy (2.7) and (3.2) for any v€ X N ChJ( 0,1). Using these 
equations it is easy to see that 

<uj(t).v> + <q N Du N (t),Dv> - <F N (t),v> - (<^ P N u(t),v> + <q N DP N u(t),Dv>) 

= 0 - (<^r P N u(t),v> + <q N DP N u(t),Dv>) 

+ <u t (t),v> + <qDu (t) ,Dv> - <F(t),v> 

where F(t) = F(t;r), F (t) e F(t;r^), and q^ = ^ + H • It thus follows 
that 

< (u N (t) - P N u(t)),v> = - <q N D(u N (t) - P N u(t)),Dv> 

+ < -fit (u(t) - P N u(t)),v> + <qDu(t) - qVu(t),Dv> + <F N (t) - F(t). 


v>. 
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Letting v = u"{t) - P u(t)€ X^, we argue that 


|u N (t)-P N u(t)| 2 < - ra|C(u N (t)-P N u(t))| 2 

+ (u(t)-P N u(t))| |u N (t)-P N u(t)| 

+ ((2m) - ** |qDu(t) - q N D(P N u(t))|)((2m) % |D(u N (t) - P N u(t))|) 
+ |F N (t) - F(t) | |u N (t) - P N u(t)| 

(u(t)-P N u(t))| 2 |qDu(t)-q N D(P N u(t))] 2 

|F N (t) - F(t)| 2 t |u N (t) - P N u(t)| 2 , 

12 2 

where we have repeatedly used the inequality ab <_ (a + b ). Defining 
w N (t) s |u N (t)-P N u(t)| 2 , the above estimates reduce to 

w N (t) - 2w N (t) < | £ (u(t)-P N u(t))| 2 + i |qOu(t)-q N D(P N u(t))| Z 

+ 1 F N ( t ) - F ( t ) | 2 , 


so that an application of the Gronwall inequality yields 


w N (t) 


i N ... d N ,, u 2 2T , N. N , N. N, 

|u (t) - P U ( t ) | _< e {-T 1 + T£ + T3 + T4} 


where 


N _ I ,,N/ n Nx d N , n -x ,2 
t-| - I u ( 0, y ) - P u(0;y) I 

t 2 = (u(s;y) - P N u(s;y))| 2 ds 

T3 = / | qDu(s,y ) - q N D(P N u(s;Y))| 2 ds 

Tw = / |F(s; r N ) - F(s;r)| 2 ds. 

4 0 
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N 

From hypothesis (H3), it follows that + 0 as N + ». We are also able to 
argue the convergence of to 0 since t 1 ^ = |P N Uq - P N u 0 1 £ |ujj - u Q |. 

To consider t 2 , it is useful to note that, for v fixed in X^, the function 
defined by g^(t) e <u(t) - P^u(.t),v> is identically zero (from the definition 
of P N ) so that 0 = gj(t) = <u^( t) - P N u(t),v>. But this is true for every 
v £ and all t € (0,T), so it must follow that — P^u( t) = P^u t (t). Thus, 

x p = / IM S ) " P N u t (s)| 2 ds 

where, for each s e (0,T), the integrand converges to zero as N + this claim 
may be verified using (3.4) and the fact that u^.(s) € L 2 (0,1) and 
V- is dense in l_ 2 (0,l). The integrand is dominated by 2 [ u t ( s ) [ 2 where s -* u t (s) 
is in L 2 ((0,T), L 2 (0,1)) (Theorem 2.1). It follows that t 2 -* 0. 

Finally, for N sufficiently large, 

mr” < / |(q-q N )Bu(s)| 2 + / T |q N D(u(s) - P N u(s))| 2 
J 0 0 

T T T 

i 2 / |(^-- ^)Du(s)| 2 + 4 / |H- (cL - 4 > 2 )Du(s) | 2 + 4 / | (H- - H J<}> 2 Du(s)| 2 
0 1 1 0 ^ ^ 0 5 5 

+ rfi 2 / |D(u(s) - P N u(s ) ) | 2 
0 

2 2 T 2 J 

+ I *?" #2 I J / l»u(s ) 1 2 + 2 m \l - ? N | / |Du(s) 1 2 

+ (mc 2 N -1 ) 2 / |A(q)u(s)| 2 
0 

T o 

where we have used (3.5) in the last inequality. Further, / |Du(s)| = 

0 

T -1 T - 1 T 

/ <Du(s), Du(s)> £ m / <qDu(s) ,Du(s)> £ m / |A(q)u(s)| |u(s)| < » 

0 0 Q 
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(since Ug € V-) so that + 0 as N -* «. 

Finally, it is possible to lift the requirement that Ug belongs to V- and 
to prove a more general convergence theorem. 

Theorem 3.2 . Assume {y N } is given in r, y N = (^.^^r^Ug) , such that 
y N -»■ y = (i,^i , 4 > 2 ,r,Ug) £ r. Then, for each t e (0,T), 

u N (.t;y N ) - u(t;y) in L 2 (0,1) ' 
as N -> “. 

Proof . Let e > 0 be given and define q = ^ + H-^» = and 

5 = ( S »^i »^2 *^) * Si pce V q 1S dense in Lg (0 , 1 ) , there exists ip e V- sufficiently 
close to iig € L 2 ( 0,1 ) so that we may argue that 

|u N (t;y N ) - u(t;y)| £ |u N (t; (s N ,Ug) ) - u N ( t ; ( s N ,ip ) ) | 

+ | u N (t; (s N ,i/>) ) - u(t;(i,^))| 

+ |u(t;(s,i/0) - u(t; (5,Ug) ) | 

< e » 

for N sufficiently large. 

N 

Here we have used the continuous dependence of u , u on initial conditions 
(uniform in s^, i, and N), the inequality |ug - <_ |Ug - Ug| + |ug - ip | , and 

the findings in Lemma 3.2. The proof of the theorem thus obtains. 
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To this point, we have focused on state variable convergence (of to u) 

once the convergence of any sequence of parameters has been guaranteed. In 

reality, however, we have yet to establish whether any sequence of solutions 
-N 

(y } for the approximating estimation problems is indeed convergent; even then 
we have no assurance that the limiting parameter y is in fact a solution to the 
original parameter identification problem. In our next theorem we consider this 
problem and indicate when an approximate estimation problem may be used numer- 
ically to compute an approximate solution for the original problem. The proof 
of the theorem is similar to ideas found in [12], [14]. 

Theorem 3.3 . For each N let y^ denote a solution for the problem of minimizing 

N * Nl m 

J over r. There exists y € r and a subsequence {y K } of {y } such that 

. -Ni. * _ 

(i) y -> y in the product topology on r, 

(ii) u^ k (t;y^ k ) -* u(t;y ) for each t£ (0,T), 

/ x Ni/ _Nk * 

(iii) J (y ) -*■ J(y ), and, 

(iv) y is a solution to the original parameter estimation problem, 
namely that of minimizing J(y) over r. 

Proof • Parts (i)- (iii) are immediate consequences of hypothesis (HI), 

and Theorem 3.2. To prove part (iv) it suffices to note that 

, *, Nk -N|e 

J(y ) = lim J k (y k ) 

N^-*- 00 

Nk 

£ lim J (y) 

N r~ 


= J(y) 
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-Nk Nk * 

for any y € r (y is a minimizer for J over r), so that y is a solution for 

the problem of minimizing J over r. 


Remark 3.1 . Although it has been assumed in our discussion so far that p = 1, 
there are no difficulties associated with extending our ideas to more general p 
(p = k-j + H^.^, (k.|,k 2 ) G K = ( (k-| jk^) £ C[0,1]xC[ 0,1] | 0 < m <_ _< m}) . 

For example, it is easy to see that the arguments used in the proof of Theorem 2.1 
change very little if, instead of A(q) and the l _2 inner product, one considers 
the operator B(p,q) = A(q) and the weighted inner product defined by 

<4>t ?>p = / p\K- For each N and y^ = given in r (r a 

compact subset of $ x LgCOJ ) x k), we may define approximating equations in the 
variable u N (t) G X N (q N ) by 


<p N u^(t),v> = - <q^Du^(t) ,Dv> + <F(t;r^),v> , veX^(q^) 

N, n x d N N 
u (0) = P u Q , 


N N N N N M N 

where p = k-j + H ^ anc * q = <j>-| + H To establish convergence of u to u 

N ^ £ 

(as y y) one may make a simple modification in the proof of Lemma 3.2 to 


argue that 


m|u N (t) - P N u(t)| 2 < |(p N ) % (u N (t) - P N u(t) ) | 2 


< ce 2T (x7 + ^ + 4 + |5u t -p N P N u t | 2 ) 


N 

where the t. are unchanged from that lemma and the last term converges to zero. 
Because all other estimates remain unchanged, we are able to derive an analog 
to Theorem 3.3, i.e., we are able to treat the case of a general (possibly 
unknown) parameter p. 
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Remark 3.2 . To this point, we have developed an estimation theory based on state 
variable approximation only; that is, we used u N to construct an approximate 
fit-to-data functional J which we then tacitly assumed could be minimized 
(numerically) to obtain an "optimal" y N €, r. Of course, one cannot actually 
use a computer to implement such a parameter search since r is in fact a 
functional parameter set (y £ r contains the functional components , 4 > 2 » r ,UQ, 

K-|,and < 2 ). The problem of further approximating the parameter set r has been 
the subject of recent studies (see [10], [14], and [23]); we shall summarize, 
in particular, the results of [14] as they pertain to the problem at hand. 

For ease of presentation we shall assume that r, ug> > and k ^ are known 
(there is an easy extension of these ideas to the case where these functional 
parameters are unknown) so that y = ( 5 , <j>-|, cj^) is the vector of parameters to 
be estimated. Since we use cubic B-spline approximations to approximate the 
functional parameters in our numerical examples (Section 4), we shall restrict 
our attention to a theory based on cubic splines only; a more general theory 
may be found in [14]. To this end, we take the (more regular) parameter set r 
to be a subset of 

S(m,ifi) = {s = (c, 4>-j , 4 > 2 ) € [ 6 , 1-6] x C 1 [0, 1 ] * 

[ 0, 1 ] 1 0 < m £ (x) <_ m for x £ [ 0 , 1 ], 

4> i e H 2 (0, 1 ) , and Id 2 ^! < m , i=l, 2}, 

(6 £,(0,1) is fixed) and assume thatris compact in the R* c"*[0,l] xC*[0,l] topology. 
For each M we define the finite-dimensional (approximate) parameter sets 
I'M by = i^(r); here i^ : R k c"*[0,1] * C^[0,1] — > R x C 2 [0,1] x C 2 [0,1] 



26 


is given by i M U, f,) e ( 5, x M,? ^ , l M ’ ?! f> 2 ) where I M ’S is defined to 

be the (unique) cubic spline function i>(x) satisfying 'Hxp = ^(x^.) » 

k = 0,1, .... 2M, and Df(xg) = D<f>(xg) , D<f>(x 2M ) = (see, for example, 

Chapter 4 of [29]). The knots xJJ are the 5 - dependent knots described earlier 

M M 

in this section, i.e., x^ = k£/M, k=0, ..., M, x^ = E, + (k-M)(l-s)/M for 

k = M+l , . . . , 2M. We remark that we are also able to construct a parameter 

approximation scheme based on a uniform mesh (of mesh length 1/M) for components 
MM M M M 

<f>-|, <f> 2 of (5, 4>-| , <j> 2 )6r ; however, as is true with the approximation of state 
variables, the resulting numerical scheme is greatly simplified if the mesh 
depends on 5, as well as on M. We shall defer to Section 4 a more detailed 
discussion on computational features of the resulting algorithm. 

It is not difficult to use the ideas of [29; Chapter 4] to argue that, for 
fixed M, the mapping (c,<t>) -> I M ’ 5 <}> : [6, 1-6] x ^[0,1] -> C[0,1] is continuous 
and thus r is compact in the R * C[0,1] x C[0,1] topology. In addition, given 
e & (0,m), we may use a variation of [29; Theorem 4.5] to see that, for M > m 
sufficiently large (the choice of ® is independent of 5), 

Me 

0 < m-e <_ I ’ ^ <j> <_ m + e , 

2 2 - 

for all 5 €> [6, 1-6] and all <j> e H (0,1 ), |D <j> | <_m, m<_ <f> (x) <_m. Therefore, for 

M 

M sufficiently large, r is a parameter set satisfying all conditions needed 

(namely r c 5 (m-e, m+e) and r satisfies hypothesis (HI)) order to apply 

the parameter estimation/state approximation theory developed thus far. In 

-N,M 

particular, for each N and M _> til, there exists a solution y to the, problem 
of minimizing over r M . From the construction of r M (and the compactness 
of r) we know that there exists a sequence {y ’ } in r with y ’ = i (y ’ ) 
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and a subsequence {y 


N k ,M J 


} such that y 


A.* 1 ' 


■> y e r in the R x C[0,1] x 


C[0,1] topology. In addition, y^’ M satisfies 


jN k( -N k ,MJ) < jNk( Y)> y6r f 


from which it follows immediately that 

(3.8) jNk(-Nk.M J ) i jNk </(})) , . 


Applying arguments very similar to those in [14], we may use the convergence 
.M 

of |i (y) - y — > 0 as M — > oo s uniform in y & r [29; Theorem 4.5] to see 
that |y^’^ “ Y -> 0 as N^M^' -> <» and that, passing to the limit in (3.8), 

(3.9) J( Y *) < J(y), y e r. 

★ 

The parameter y is thus a minimizer for J over r. We summarize these findings 
below. 


MM - N M 

Theorem 3.4 Let r = i (r) and let y ’ denote a solution to the problem of 


. . . N M 

minimizing J over r . Then there is a subsequence {y 


N of f^.M } such 


that Y Nk ’ MJ -> 


Y Pk " — > y » where y is a solution to the problem of minimizing J over r. 

In fact, any convergent subsequence has as its limit a solution to the original 
estimation problem. 


3.1. Approximate Estimation Problems Associated with "Discrete" Data 

It is possible, under additional smoothness assumptions on solutions, to 
use variati onal-type estimates (similar to those found above or in [10], [30]) 
to argue pointwise (in x) convergence of state variables; i.e., u^(t,x;Y^*) 
u(t,x;y) whenever y^ + y» for (t,x) G (0,T) x [0,1]. Results of this type 
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lead naturally to a statement about the approximation of a solution for the 
problem of minimizing the "pointwise" fit-to-data criterion J (see (2.3)) 
over r. We shall briefly summarize our findings below. 

It is not surprising that for this formulation we require additional 
smoothness assumptions on the parameters. Specifically we take y = (s,Uq) £ r 
where r is a subset of r = 3 * Hq( 0,1 ) , and 3 = (s = (g,^ ,<j>2»r) € 3 | 

|J ) 4> i L 1 1 = 1 >2}. 

In addition to hypotheses (HI) - (H3), we make the following assumption 
(which in general may impose additional constraints on parameters and the 
applied force f): 

(H4) For any y £ r, the mapping s -> u^.(s;y) is in l_2((0,T); hJ(0,1)). 

N/ x 

Defining approximate state spaces X (q) as before, we seek approximations 
U N to u, where, for any given y = (s,^ ,<f> 2 , r, u Q ) e r, q = 4> 1 + H * 2 , 

u^(t) = u^(t;y) satisfies 


r <u^(t),v> = - <qDu N (t),Dv> + <F(t;r),v>, te (0,T), 


(3.10) { 


N 


u"(0) = P N (q)u 


for all v€ X N (q); P N differs from defined in (3.2) in that P N : Hq ( 0 , 1 ) •> X N (q) 
is the orthogonal projection in the hJ( 0,1) (rather than L 2 (0,1)) topology. 

It is not difficult to see that there exists a unique solution u N (y) of (3.10). 

In addition, for each (t,x) £ [0,T) x [0,1], the mapping y -> u^(t;y)(x) is 
continuous (in the 3 x Hq( 0,1) topology on y). 

The convergence result that follows is a pointwise analog of Lemma 3.2. 
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Lemma 3.3 . Let {y^} be given in r such that y^ y in the 3 x Hq ( 0,1) topology, 
Y = (I,^i r ’ ^ ^0 ^ ^q» t * 1en u ^(t;y^) -*■ u(t;y) (as N -»■ «) in the 

Hq ( 0,1) topology, uniform in t 6 [0,T); here u^(y^) and u(y) are solutions of 
(3.10) and (2.7 ) corresponding to y^ and y, respectively. 

Proof: We shall write u N = u N (t;y N ), u = u(t;y), q N = + H 

q = ^1 + H-^ 2 > = F(t;r^), and F = F(t;r) throughout. We note that 

|D(u N - u) | £ |D(u N - X N u)| + |d(I N u - u) | 

so that, using (3.6) and the fact that u(t) € V- for t £ 0, it suffices to 
show |d(u N - I N u)| -> 0 as N -* ~. 

Using (3.10), (2.7), and arguments similar to those in the proof of 
Lemma 3.2, we may argue that, for v € X N = X N (q N ), 


N„N NNNN N - NN 

<u t - I u t ,v> + <q Du - q DI u,Dv> = <u t - 1 u^, v> + <qDu - q DI u,Dv> 


, C N c 
+ <F - F,v> 


and in particular, using v = u^ - I^u^., 

luj - I N u t | 2 + <q N D(u N - I N u), D(u^ - I N u t )> 
£ J |u t - X N u t 1 2 + |uj - I N u t 1 2 


+ \ | F N - F I 2 + <qDu - q N DI N u, D(ujj - I N u t )>. 


We thus find that 
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1 d 


j -£ {|(q'T 2 D(u N - X N u) | 2 - 2<qDu - qVu, D(u m - !%)>} 


N..N, 


.N -N. 


- \ l U t ' lNu t ! 2 + \ l pN - F ! 2 - <ql>u t - q N DI N u t , D(u N - X N u)> 


d N N 

where we have used the fact that ^ 1 u = I u^.; it therefore follows that 


m|D(u N - I N u) (t) | 2 < ^(t) + ^ a^t) + a^O) + (ffn-lh^O) + a 3 


(3.11) 


N N n 
+ a 4 + a 5 + a g , 


where 

a"(t) = |q»u(t) - q N DX N u(t )| 2 , 

a”(t) = |D(u N - I N u)(t )| 2 , 

N _ r T | 1 2 

a 3 ^ I U t “ 1 U t I ’ 

N _ r T I j-N .-1 2 

a 4 / l F “ F l 
0 

<75 = - 2 / <(q - q N )Du t , D(u N - I N u)> , 

3 0 Z 

and 

a g = _ 2 J T <q N D(u - lV), D(u N - I N u)> . 
o 0 c c 

Using (3.6) and arguments like those for in the proof of Lemma 3.2, we 


find that 
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a^t) < 2|(q - q N )Du(t)| 2 + 2m 2 |p(u - X N u) (t) | 2 
1 2k(N) | Du ( t ) | 2 + 2m 2 (nmN)‘ 2 i A(q)u(t) | 2 , 

where 

k(N) = 4(|^ - cf^ 2 + |^ 2 - 1 2 ) + 2m |£- £ N | -> 0 

Thus, for t £ [0,T) , 

°?(t) 1 2k(N) sup | Du(t ) | 2 + 2m 2 (ni7rN)" 2 sup 

t<0,T) t<0,T) 


Considering aJj and <jg, we find 

a!j(t) ^ sup | D(u N - X N u) (t) | 2 
te[0,T) 

and 


a 


N 

5 


0 z 
i^k(N) / 1 Du t | 2 + f 


♦ * /W - tW 

0 

sup |d(u N - l N u) ( t ) | 2 
te[0,T) 


In addition, we may integrate by parts to show that 

.N 


2N T x 


>6 - - 2 Z / / k »(u t - I N u t )q N P(u N - X N 
k=l ON r t 


u) 


'k-l 


2N T x? 

= 2 Z / / k K 

k=l ON r 
x k-l 


d S N -> 00. 


| A(q)u(t) I 2 


X N u t ){Dq N D(u N - I N u) + q N • 0 } 
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< m^mTrrf 1 /V - X N u t | 2 + gjjfr J T |»(u N - I N u)| 2 } 

< 8m 2 Tm _1 / V - I N u.| 2 + § sup |D(u N - I N u)(t) | 2 , 

0 L Z a u[o,T) 

N 

where we have used |D<ji. | < m, i = 1,2. Therefore, from (3.11), 

I sup |d(u N - I N u)(t)| 2 < (2m -1 + l)a!J + (rm-T ) |X)(P N Uq - X N G n ) | 2 
8 te[0,T) 1 00 

+ (8m Tm ^ + 1 + cr^J 

1 T 2 

+ 4Tk(N)m" / | Du . | L , 

0 z 

where a^, aj + 0 as N ■>• «. We may apply (H4) and standard spline estimates 

(see, for example. Theorem 2.4 of [29]) to also show that ■> 0 as N -> °°. 

It remains to consider |d(P^Uq - I^Uq)| 2 : We note, using the properties of 

N 

P , that 

|d(P N Uq - I N u q )| < |d(p N Uq - P N u 0 )|+ |d(p n u 0 - u 0 )| + |l>(u 0 - l N u 0 ) | 

< | D(Uq - u 0 )| + 2 |d(1 N u q - u 0 )|, 

N 

where each term converges to 0 as N from the convergence of Uq to Uq in 
H J and (3.6). The proof of the lemma thus obtains. 


Finally, we prove a more general convergence result. 
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Theorem 3.5 . Let {y N } be given in r with y N y e r in the $ * hJ(0,1) topology. 
Then, for each t € [0,T), 

u N (t;y N ) - u(t;y) in hJ(0,1) 

as N -*- <». 

Proof : We first demonstrate the continuity of the mapping 

u Q + u(t; (s,Uq)) : hJ -> hJ, uniform in te [0,T) and $. Let u Q , u^ £ hJ 

and u = u(t; (s,Ug)) , u 1 = u(t;(s,Ug)). For any v£ hJ, 

I 

<u t " u t’ v> = " <c l( Du " Du'),Dv> 
so that, using v = u t - uj. € hJ, 

|u t - \\ 2 = - 7 ^ |q % o(u - u' ) I 2 . 

Applying the Gronwall inequality, we find that 
|q % D(u - u')(t)| 2 < |q % D(u 0 - uj)| 2 , 

or that 

| D(u - u * ) (t) | 2 < mm _1 |D(u 0 - Ug)| 2 , 
so that the continuous dependence result obtains. 

We may construct similar arguments to demonstrate (using (3.10)) that 
|D(u N (t;u Q ) - u N (t,Ug))| 2 < mm" 1 |d(P N ii 0 - P N Ug )| 2 

<_ mm -1 |d(u q - U g)| 2 
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so that the mappings u Q + u N (t;(s,u Q )) : hJ -»■ hJ are also continuous, uniform 
in t, s, and N. 

Let {y^} be given in r with y^ y, where y = ,if>2,r,UQ) G r. To 

argue the convergence of u^(t;y^*) -> u(t;y) in Hq for arbitrary Uq G Hq (using 

estimates like those in the proof of Theorem 3.2) we need only demonstrate that 

V-, q = ^ + H-$ 2 , is dense in hJ. To this end, let e > 0 and ip be arbitrary 

in Hq ( 0 , 1 ) and define ip = ip - h where h is given by 


h(x) 


ip(c)x/?, 


4 


P(x) 


x G [ 0 , 5 ] 
x € (?,!] . 


Here p is a quadratic polynomial satisfying p ( 1 ) = 0, p(5) = ip(i)> and 
Dp(s) = q(I“)ip(I)/iq(I + ) (it is easy to show such a p exists, under the condition 
5 / 1 ). It is clear from the construction of h that we have hG V-, ip G Hq(0,1), 
and i(x) = 0 for x = 0,5,1, so that ip G hJ(0,s), ip G hJ(c,1). From the 

definition of hJ, there exists ^ G Cq(0,5) and ^ G Cg(5,l) such that 
2 ~ 

U - Cl I < e/2 and |ip - 5 J 2 < e/2. Finally, defining 5 on [0,1] 

Hq(0,5) 2 hJ(c,1) 

by 5 = + h on [ 0 , 5 ], 5 = ? 2 + h on (1,1], it is easy to see that 5 G V- 


(5 G Hq ( 0 , 1 ) , qD 5 G H^ (0,5) , qD? G H 1 (1, 1 ) , and qDs is continuous at x = 5 ) 


and |ip - ? 1 


Hq ( 0 , 1 ) 


= k 


hJ(0,5) 


lip - 5 0 r, < e. Therefore, V- is 

2 Hq( 5,1 ) q 


dense in Hq(0,1) and the proof of the theorem is complete. 

Finally, we return to the problem of (approximately) determining a 
minimizer y* for the (pointwise) least squares criterion J. The proof of our 
final result uses the estimates derived in this section, following the proof 


of Theorem 3.3. 
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Theorem 3.6 . For each N, let denote a solution for the problem of 


^ N ~ m n n »i 

minimizing J over r, where iT(y) = l l |C(t.,x.;y)u N (t.,x.; 

i=l j=l J 1 J 


iv) - u.. 


U 




and u is the solution of (3.10) associated with yer. Then there exists y e r 

-Nk ~N ~Nk ~^k ~Nk ~ 

and a subsequence {y } of {y } such that y y , J (y ) ■> J(y ), and y 

is a solution for the problem of minimizing J over r. 


One may also easily modify the arguments given above to include p as a 
parameter (see Remark 3.1) and to prove a "double approximation" result 
similar to Theorem 3.4 for both the state and estimated (functional) parameters. 
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4. Implementation and Numerical Findings 

A desirable feature of the spline-based scheme developed in preceding 

sections is the ease of implementation of the approximation ideas, especially 

when the points of discontinuity ^ , i=l,...,y-l, for coefficients are unknown 

and to be estimated. In what follows we describe how the particular state 

approximation framework chosen here serves to facilitate (from a computational 

standpoint) the parameter estimation/approximation process. We conclude the 

section by presenting our findings for some representative test examples. 

We begin by examining the approximating ordinary differential equation 

(3.1) rewritten here in terms of w^(t;y) = (w^(t;y), w^tiy), . W 2N-1 ( t;y ) )"*"> 

N 

where the w^, defined in Section 3, are the coefficients in the expansion 
N 2N-1 N N 

u w (t;y) = l w^(t;y)B- (q). Using this notation, the ODE may be written 
i=l 


(4.1) 


'W(t) = - K N w N (t) + G N (t) , tfc(O.T) , 

^w N (0) = Wq ; 


here the (2N-1 )-square matrices Q 


Q N (y) and K N = K N (y) have entries 



<Bj (q), B^(q)> , 


K i,j = DB i( c i) > » 

while the perturbation term and initial condition satisfy 

G N (t) = G N (t;y) = (<F(t;r), B^qK ..., <F(t;r), (q)>) T 

and 

w 0 = w 0 (y) E "*( <u o* B^*(q)>, ..., <u Q , B^ N _i(q) > ) > 


respecti vely. 
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Implementation of a computational scheme to estimate (approximate) para- 
meters in (4.1) is greatly facilitated by the choice of basis elements for the 
N approximate state space X N (q). In order to best indicate some of the advantages 
of this approximation framework, we first consider the special case where only 
q is unknown, where q = <j>-| + H^, and and <j> 2 are constants. First, it is 

easy to see how our choice of a linear spline approximation scheme yields matrices 
N N 

K and Q that are quite simple in structure: For a given value of q, the inner 

products appearing in these matrices may be determined from explicit formulas 
(depending on N and £), a few of which are given here. For example, diagonal 
entries in the (tri diagonal ) matrices and are given by 

Q-j = 2S/3N, i = 1, ..., N-l , 

Q N % = 1/ 3 N , 

Q i,i = 20-S)/3N, i = N+l , ..., 2N-1 , 

(4.2) K^. = 2^/5, i = 1, .... N-l , 

(4 - 3 ) k n, n = N( h /5 + N(< h + + 2 )/(i" 5) ’ 

and 

(4.4) K^. = 2N ( 4> 1 + * 2 )/(l-e), i = N+l, ..., 2N-1 , 

with similar representations for off-diagonal elements. We note that we are 
able to avoid time-consuming and error-producing numerical quadratures; in 
addition, our approach is more desirable (from a computational point of view) 
than a method based on a uniform mesh size. For example, if for each N we simply 
subdivide [0,1] into units of length (so that position of e , is not taken into 
account) the matrix Q will be fixed throughout the estimation process; this 
however is at the expense of considerable added difficulties associated with 
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N 

evaluating entries in K . Using a uniform mesh, some of the inner products 
must be "broken up" at the point s, e.g.. 


,N 


,N 


N _ n N 


<qoBj , db?> = 4> 1 J 0 DB‘: pb': + (* 1 + * 2 ) 


1 


db'J dbJ 1 

J * 


requiring (multiple) numerical quadratures every time that q (and thus 5 ) is 

updated. In contrast, with the ^-dependent structure chosen here we need only 

recombine simple algebraic expressions (such as those given in (4.2) - (4.4)) 

N 

to obtain the elements of K . 

Many of these computational advantages are still present in the case where 

the <f>. are not assumed to be constant. If, for example, M and N are fixed 
M 

and r consists of cubic spline element approximations for <f>^, ^ (defined 

M 

^-dependent mesh of points x^ ; see Section 3), many of the quadratures may 

still be performed in advance of the iterative process. In particular, if we 
• M/..X M Ji /..x o ........ J1 


let <j> n (x) = 


m=l 


n,m m 


G (x) , for n=l,2, where G m are the usual cubic B-spline 


m 


M 


basis elements defined using the mesh points (x k , k=0, ..., 2M} , we find that 


N M N N 

K. • = <q DB., DB.> may now be written as 
1 >J J 1 

k(M) 


k(M) 


(A n _ V- M J1 n „N n oN x v / M . M v J1- D N . d N 

(4 ’ 5) K U " Y l,m <G m DB j’ Jn (y 1 ,m + Y 2,m } <G m DB j ’ DB i > L 2 (e,l )’ 


Since simple explicit algebraic expressions (in terms of £ and M) exist for 
M 

G m , the quadratures in (4.5) may also be worked out easily in advance (yielding 

expressions involving 5 , M, and N); as q-iterates are updated (i.e., coefficients 

M » N 

y are updated) it becomes a simple task to calculate the new entries in K . 
n,m r r 

★ 

We consider here numerical examples where y is known and we have generated 
synthetic data for use in testing our ideas. In all examples presented here, 
we assume that r and Uq are known and fixed at their true values so that only 
q = <f>i + H^ 4>2 is unknown (i.e., y = ( 5 , <j)-|, (j^)) and to be determined. The special 



39 


problems associated with estimating this discontinuous coefficient have been the 

focus of our efforts throughout; the problem of identifying continuous functional 

parameters and initial conditions has been considered elsewhere [13], [14], [16]. 

For each example that follows, both y and u(y ) are selected in advance while 

the appropriate forcing function f is artificially determined by substituting 
* * 

y , u(y ) into (2.1). For chosen sample times t^ , i =1 , ..., n, and sampling 
locations x. , j=l, ..., n (discrete data is used for these examples), data 
is generated by setting u^- = u ( t^ , x^.; y ) , with random noise added in some 
cases. We note that the sample data is not generated using our 
spline-based scheme; rather, the data is constructed from an analytic expression 
for the solution and thus is independent of the methods we illustrate here. 

We begin the parameter estimation process by supplying an initial guess of 

0 

Y to IMSL's minimization routine ZXSSQ (a Levenberg-Marquardt algorithm) which 
numerically attempts to determine a minimum, for given N, to (using e = 1 
in (2.3)) over a fixed constraint set r M . Here u N (y) is the solution to (3.1) 
calculated using IMSL's DGEAR, an ODE solver, where the known values of Uq and 
f are used in the equations. We note that although we are actually using the 
cost functional associated with discrete observations u-., the approximating 
equations (4.1) differ somewhat from those defined in Section 3.1: Indeed it 

is not surprising that, in practice, we obtain pointwise convergence of the 
approximating states under hypotheses more general than those needed in Section 
3.1 so that we may, in fact, relax some of the restrictions on the approximating 
system. 

Example 4.1. In our first example we take 

* I 15. , 0 < x < .6 

q (x) = { 

I 50. , . 6 < x < 1 
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★ * 

and define u(t, •; y ) e dom A(q ) by 

* fx(70-100x)(t 2 +2) , 0 ^ x < .6 

u(t,x; y ) = ? 

L(15-15x)(r+2) , .6^x<l. 

3 

In examples 4.1 .a - 4.1 .c below we seek to estimate y = U, 4>-j , ^^rcR 

r\ 

(with true value y = (.6, 15., 50.)) using an initial guess of y u = (.8, 30., 30.). 

N 

In each case we obtain the converged values y for N=4, 8, 16, and 24, using 
y° to start the iterative scheme for N=4, and previous converged values as start- 
ups for N=8, 16, 24 (e.g., y^ is used as initial guess for the N=8 run). 

★ 

Example 4. 1 .a. Data is generated for this example using u.. = u(t. , x.; y ) 

• J * J 

for t. = . 5i , i=l, ..., 4, and x- = . 1 j , j=l, ...» 9. Our findings are reported 

* J 

in Table 4.1. a. 

Example 4.1 .b. We repeat the last example except that spatial sampling locations 

are now given by x- = .lj + .05, j= 0, 1, ..., 9 (so that there is no spatial 

J 

* 

observation point at ? , the point of discontinuity). We summarize our results 
in Table 4.1 . b . and note that there is little change between this example and 
Example 4.1 .a. 

Example 4.1 .c. We repeat Example 4.1. a, but add noise to the data. In this 

* 

case we define u. • = u(t. , x.; y ) + r . . where { r i - } are Gaussian random numbers 

which (with 98% certainty) fall in the range [-.060, .060], u = £ u-./(nn). 

i . j J 

Our findings for this example are summarized in Table 4.1.c. 

In the examples that follow we shall shorten our discussion by 
abbreviating the length (and number) of tables and by displaying some results 
graphically. The rather detailed presentation given for Example 4.1 was pro- 
vided simply for the purpose of observing if noise in the data or changes in 
the placement of data affected the outcome. 
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Table 

4.1. a. — 

Example 4.1 .a. 



N 

_N 

_N 

<h 

_N 

*2 

~N 

J 

CP time 
(secs) 

No. of 
i terates 

4 

.623 

14.669 

51.950 

1.5 x 10 2 

28. 

13 

8 

.602 

14.845 

50.672 

1.5 x io° 

54. 

7 

16 

.600 

14.961 

50.095 

8.8 x 10‘ 2 

202. 

7 

24 

.600 

15.000 

50.000 

5.6 x 10" 9 

141. 

4 



Table 4.1 . b: — 

Example 4.1.b. 



N 

-N 

5 

-N 

*1 

-N 

*2 

~N 

J 

CP time 
(secs) 

No. of 
iterates 

4 

.621 

14.956 

48.494 

9.0 x 10 1 

32. 

20 

8 

.607 

15.009 

50.063 

3.8 x 10 1 

35. 

7 

16 

.601 

14.991 

49.728 

1.2 x io _1 

355. 

13 

24 

.600 

15.000 

50.000 

5.7 x 10" 9 

239. 

5 



Table 4.1 .c. 

-- Example 4.1.c. (Noisy data) 


N 

f N 

*1 

-N 

9 2 

o N 

CP time 
(secs) 

No. of 
i terates 

4 

.621 

14.730 

51.573 

1.6 X 10 2 

27. 

10 

8 

.599 

14.887 

50.434 

8.1 x 10° 

68. 

10 

16 

.598 

14.991 

50.296 

8.0 x 10° 

178. 

5 

24 

.597 

15.006 

50.149 

8.3 x 10° 

733. 

8 
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In Example 4.2 below, we illustrate the use of our methods in problems 
with two discontinuities ^ , £ 2 , in q; the example also serves to illustrate that 
we are able to accurately estimate even when the forcing function f does not 
contain discontinuities at each of those points. 

Example 4.2. We seek here the "true" value of q given by 


* 

v 1 

O 

O 

X < 

.2 

q = 

1 6.0 , .2 < 

X < 

.6 


^0.5 , .6 < 

X £ 

1. 

In this case, y 

* * * 

* 

* 

ill 

vT"< 

vPT 

PO 

s* 

> 

<j>2 


true solution corresponding to y is 


6. , .5) and the 


. f 30 * 

u(t, x; y ) = < 5x + 5 

l-200x 2 + 300x - 100 , 


0 £ x < .2 

.2 _< x < .6 

.6 £ x £ 1 , 


with data available at t. = . 5i , i=l,..., 4, and x- = .1- , j=l,..., 9. 

I J J 

A sample of our findings is given in Table 4.2 below, where the converged values 

~8 _6 

reported were obtained after 501 CP seconds, with >3 = 8.3 x 10 . 



Table 4.2 - 

Example 4.2 




?! ? 2 

♦l 

<f>2 

*3 

Initial guess:. 

.300 .700 

5.000 

5.000 

5.000 

Converqed values 
(N = 8): 

.200 .600 

1.000 

6.000 

.5000 

True values: 

.200 .600 

1.000 

6.000 

.5000 
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-jc -fc ^ 

We consider now two examples where the "true" q = <f>i + H *$2 involves 

"k 

nonconstant values of ^ and ^ . In each case we search for approximate 
(j >1 and <j >2 in the cubic spline space constructed using an M =1 level of approx- 
mation (see Section 3). 


* 

q = 


Example 4.3. Here we seek to estimate the "true" parameter 
: + 12 , 0 £ x < .6 

J100x 2 /9 , .6 £ x < 1 

starting from the initial guess for q of q 0 ' = 3 on [0,1] (with start-up value 
for 5 of K = -5). The solution 


u(t, x; y ) 


"(70x - 100x 2 )(t 2 + 2) , 
2 

(15(l-x) (t + 2) , 


0 £ x < .6 
.6 £ x < 1 , 


is used to generate data at ^ = .5i, i=l, 4, and x^ = .1., j=l , 9. 

-N M N M N M 

In Figure la we compare the estimated q ’ = ^ + H-n (N=16, M=l) 

★ 

with the "true" coefficient q . Figure lb is the same graph that has been 
enlarged and restricted to the interval [.4, .63] in order to better distinguish 
between "true" and approximate curves. 


Example 4.4. Again we estimate a functional parameter with true representation 
given by 

’.424 - 40x , 0 £ x < .3 

J90(x-.3) 2 + 18 , .3 £ x £ 1. 

Data is generated as in Example 4.3, using instead the solution 


* 

q = 


u(t, x; Y ) = 


200xt 4 ( . 5-x) 

l 2 . 


, 0 £ x < . 3 

J 7. 143t^(l-x) , . 3 £ x £ 1 . 

We began the parameter search with the start-up guess of 



TRUE Q 

INITIAL GUESS FOR Q 

APPROXIMATE Q 


Figure lb: Example 4.3 





0 < X < .2 


45 


18. 


(j48. , .2 <_ x £ 1 

_N,M 

(see Figure 2a) and obtained the converged value of q ’ (N=24, M=l) that is 

, . .24 * 

depicted in Figure 2b. Note that i is not approaching ? ; in fact, we 

observed that the iterates for e, never changed from the initial guess of 
0 * 

5 = .2 (recall 5 = .3) throughout the iterative process. We note that the 

software package did perform fairly well in its attempt to estimate the 

-24 _24 

approximate functional shape of the parameter q (i.e., <j>-| , <j> 2 are roughly 
★ ★ 


the same as <f>^ , <f> 2 on the intervals [0., .2] and [.3, 1.]; between x= .2 

-24 

and x= .3 there is discrepancy due to the incorrect value of 5 ). The failure 

of the numerical package to adequately estimate c may be due to some well- 
known limitations of the particular optimization scheme (Levenberg- 
Marquardt) that we chose to use with our approximation ideas: It has been our 

experience that difficulties sometimes arise when this minimization scheme is 
used to estimate more than seven or eight unknown parameters. (In this example 
there are 9 unknowns -- 5 and 4 coefficients each in the cubic spline repre- 
sentations for cf>i and <j> 2 .) That the numerical package was able to estimate 
9 parameters (and, in particular, 5) in the last example may be due to the fact 
that the difference between q(5~) and q(s + ) in that example is greater, making the 
accurate placement of 5 more critical. 

We were able to overcome the difficulties we encountered in this example 
by taking the following steps: First, observing that 5 had not changed at all 

while <j) 1 and <j> 2 appeared to have converged to reasonable values, we restarted the 

-24 -24 

iterative process holding ^ e <j>.| and <j> 2 fixed while iterating only on 

K- As is seen in Figure 2c, the converged value of 5 obtained using this approach 
-24 * 

is 5 = .299 = .3 = 5 . Finally, we "re-tuned" the coefficients in the spline 


Figures 2a-2d : Example 4.4 




Figure 2a 


Figure 2b ^ 

O' 



Figure 2c 


Figure 2d 
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expansions for ^ and ^ by iterating once again on those coefficients, this 

-24 

time holding z fixed at the new value of Z and using as start-up values for 

<J> 1 » i >2 the converged values obtained from the first iterative process (i.e., 

-24 -24 

4>-| » $2 )• The result is shown in Figure 2d. This somewhat adaptive algorithm 

to estimate accurately all unknown parameters is a common approach taken (often 
of necessity) when real data is used in connection with model -building applications 
(see, for example, [11]). 

Finally, we remark that a drawback of our approximation framework is that 
we must specify the numbe r of discontinuities in advance of the estimation 
process. Fortunately, it is possible to overestimate and underestimate this 
number and still obtain useful information. This will be the focus of our last 
two examples. 

Example 4.5. We repeat Example 4.2 except that we assume throughout that q 

is discontinuous at only one point (while two discontinuities are actually 
★ 

present in q ); we also allow spatial variation in <j>-j and ^ and approximate 

-N,M 

using cubic splines. An initial guess for q and a converged estimate q 
( N = 24 , M=l) are depicted in Figure 3 where it is interesting to note that 
the initial guess of .4 converges to a value close to that of the true 
(second) discontinuity, ^ = .6 . In addition, to the right of this point the 

• it 

estimated shape of q begins to approximate the constant function while to 
the left of that point the rapidly increasing estimated shape gives an indication 
that we have underestimated the number of discontinuities present. 

Example 4.6. We repeat Example 4.1, except that now we overestimate the number 

of discontinuities in q. We assume throughout that q = + H-. <J> 9 + H <}>.■, where 

^ 


I I i 

.a 0.3 0.4 


T 


1 

1 


i 1 r 

0.5 0.6 0.7 0.8 0.9 



true q 

initial guess 
approximate q 



Figure 3: Example 4.5 


<t> 1 , .<f> 2 » and 4>^ are constants. For an initial guess of 



r 25. , 

O 

|A 

X 

A 

on 


q° = < 

5. , 

.20. , 

.5 £ x < .7 
.7 <. x <_ 1 

9 

we obtained (N 

= 8 , 291 

CP seconds) 


8 

"14.95 , 

0. £ x 

< .503 

q = < 

14.99 , 

.503 £ x 

< .600 


.50 . 05 , 

.600 < x 

< 1 ; 


repeating the same example but with a different initial guess. 
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0 


q 


''.001 

<•001 


L- ooi 


0 <_ x < . 333 
.333 £ x < .667 
.667 < x < 1 


we observed the following converged values (N=16) 


-16 

q 


^ 2.44 , 

15.05 , 


149.88 , 


0 < x < .0001 
.0001 £ x < .6001 
.6001 < x < 1 . 


A close inspection of either result reveals that we were, in fact, able 

★ 

to accurately estimate q (as defined in Example 4.1), even though a two- 
discontinuity approximation structure was incorrectly used throughout. 


Remark 4.1 We note that all examples presented here involve polynomial or 
piecewise-polynomial state/parameter functions; from this one might conclude 
that such a polynomial structure is needed in order to effectively apply our 
spline-based methods. In fact this is not the case, as we have seen in a 
number of test examples (see [4], [7], [26] for a number of examples in the 
context of several applications). 
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5. Concluding Remarks 

In the above presentation we have given a convergence theory for algorithms 
for the special problem of estimating discontinuous functional coefficients in 
parabolic systems. We are currently working to further develop and refine these 
ideas and to extend the theory to other applications, e.g. hyperbolic (seismic) 
equations and higher order (elastic beam) systems. In particular, we are studying 
an approximation framework that imposes the continuity condition (2.5) on approx- 
imate solutions u N as well as on the original solution u. Our investigations 
in this direction involve making further (parameter dependent) modifications of 
the spline- based basis elements described in Section 3; we are also studying a 
completely different approach that involves the "tau-Legendre" ideas that have 
been successfully applied in [8] to (discontinuous coefficient) hyperbolic systems. 

We are also working to develop a related theory for two-dimensional domains, 
although for obvious reasons this is not simply a trivial extension of the ideas 
presented thus far. 

Finally, we note that we have not addressed here the problem of "identifiability" 
of parameters, or the ill-posed nature of parameter estimation or "inverse" 
problems in general. These important and difficult questions arise in all 
parameter estimation problems and are not special difficulties associated with 
the particular problem under consideration here. The reader is referred to 
Remarks 3.3 and 3.4 of [14] for further comments regarding the nature of this 
ill-posedness, of nonuniqueness and the importance of initial guesses for parameters, 

and the general unavailability of convergence rates for approximate parameters 

_N,M 
{q }. 
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